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ABSTRACT 

The Geophysical Fluid Dynamics Laboratory Modular Ocean Model (GFDL 
MOM) is used to investigate the model difference between compatible and incompatible 
surface wind and buoyancy forcing. The atmosphere is a physical system in which 
surface wind and temperature fields are related, however in most ocean numerical 
models, the wind stress and buoyancy forcing are usually specified separately, i.e., no 
constraint between the surface wind stress and surface air temperature is considered. In 
reality, only one of these two fields can be prescribed in the atmosphere-driven ocean 
models. When the surface wind field is prescribed, the surface air temperature should 
be derived, and vice versa. If the two related fields are treated as totally independent in 
forcing the ocean models the results will be distorted. Since the model solutions depend 
upon the atmospheric forcing, it is important that we study the compatibility between the 
wind and buoyancy forcings and the effect which incompatibility might have on the ocean 
numerical models. 

This study shows that the surface wind and buoyancy forcing widely used in 
ocean numerical models are incompatible. Such an incompatibility results in 21% error 
in the total northward transport of heat, 16% error in the total northward transport of 


salt, 25% error in v velocity, and 16% error in w velocity. 
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I. INTRODUCTION 


In the atmosphere the surface wind and the sea surface air temperature (SSAT) are 
closely related. As such, thermal field and the surface wind field cannot be prescribed 
independently in a numerical model. 

However, in current ocean numerical models, the surface wind and the net surface 
heat flux are usually prescribed independently. The assumption that the prescribed 
surface wind field always matches the surface thermal field (or the net surface heat flux) 
is rather unrealistic. Contrary to the assumption, the two fields do not match each other 
in many cases. The imbalance between the surface wind and the surface thermal field 
gives rise to certain adjustment processes, thus changing the original forcing field. 
Consequently, during the entire integration of the ocean numerical models these two 
forcing fields will not maintain the prescribed form. 

That there is an obvious incompatibility can be found from the pioneering work by 
Takano (1975), which states that the SSAT is a function of latitude (y) only, and that the 


surface wind stress is a function of position (x,y) and season (f) 


IA E US, Y E), VS Vx YE) (1) 


S S 


In equation (1) the atmospheric surface isotherms are given by straight lines 
running parallel to the x-axis. However, surface winds have both x and y components 


and seasonal variation. It is thus unlikely that the atmospheric surface isotherms would 


remain linear and in a steady-state when there are surface winds blowing across the 
isotherms. 

Haney (1971) states that by employing a heat budget analysis appropriate to zonally 
and time averaged conditions within the atmosphere, it can be shown that the net 


downward heat flux Q at the ocean’s surface is expressed as 
Q + O (TA = FA (2) 


where T, is an apparent atmospheric equilibrium temperature, T, is the sea surface 
temperature, and Q, is a coefficient determined from the zonally and time average data. 
The zonally averaged surface zonal winds and apparent air temperature were broadly 
used as mechanical and thermal forcing in ocean numerical models when the surface 
thermal boundary conditions (2) in Haney (1971) are met. This approach is well 
accepted by ocean modelers. 

Several combinations of surface wind stress r and T, in ocean models are listed as 
follows (Chu, 1992). In Davey (1983) model, the atmospheric temperature T, is 


specified by 


Ta SA 
Le Ue a 

Ta (y) = 7, + e cor 7 Ys <s y Ty 
Ty Vena (3) 


with 7 = 2000 km, T, = 12°C, and Ty = 8°C. The surface winds are set to zero, thus 


Cie Ve O (4) 


sf 


The important question arises as to whether the atmospheric forcing indicated by 
(3) and (4) can last 2000 days, the period over which the model was integrated. 
Geostrophic baroclinic zonal atmospheric flow is usually driven by the north-south SSAT 
gradient while the cross isothermal surface winds is generated by the resulting friction. 
Hence, the atmospheric mechanical forcing given by (4) is not consistent with the 
atmospheric thermal forcing of (3). 

There are other renowned ocean model, such as Robinson et al (1977), Semtner 
and Mintz (1977), Cox and Bryan (1984), Cox (1985), Huang (1989), and Willmott and 
Darby (1990). In general, the wind stress is purely zonal (7, = 0) and assumed as a 
sinusoidal function (or other form) of y. T, is assumed as a linear function of y, e.g., 


in Robinson el al (1977), 


Y= YN 7 
T. = T COST —————, t= 0, T, = 26 - = (y - y.) (5) 
4 0 YN Ys a i > 5 


Since the ocean numerical models are generally integrated over long time interval, 
such as 880 years in Bryan and Cox (1984), it is important to determine whether or not 
the surface wind and the SSAT as employed by the ocean numerical models are 
consistent. Where there are inconsistency, the effect of the inconsistency on the ocean 
circulation should be further investigated. 

The main purpose of this study is to address these differences. For this, 


experimentations with a three-dimensional primitive equation ocean model are performed. 


Through model simulation of the principal aspects of the Pacific Ocean, we gain a better 
understanding of the differences between the atmospheric parameters representing 
mechanical and buoyancy forcing, i.e., the surface wind and SSAT field. 

A statement of the problem in existing ocean numerical models is contained in this 
chapter. Chapter II gives a description of the Modular Ocean Model (MOM). Chapter 
II describes the compatibility between surface wind stress and buoyancy forcing. 
Chapter IV shows the errors caused by incompatible forcing. Finally, Chapter V states 


the conclusions. 


Ii. DESCRIPTION OF THE MODULAR OCEAN MODEL (MOM) 


A. INTRODUCTION 

More than two decades have passed since the first 3-dimensional primitive equation 
numerical ocean model was used to study the most basic aspects of large-scale, baroclinic 
ocean circulation (Bryan and Cox, 1967). A description of the physics and numerics 
involved was published in Bryan (1969). In this model, the prediction of currents is 
carried out using the Navier-Stokes equations with three basic assumptions. The 
Boussinesq approximation is adopted, in which density differences are neglected except 
in the buoyancy term. The hydrostatic assumption is made in which local acceleration 
and other terms of equal order are eliminated from the equation of vertical motion. 
Lastly, closure is attained by adopting the turbulent viscosity hypothesis in which stresses 
exerted by scales of motion too small to be resolved by the grid are represented as an 
enhanced molecular mixing. The temperature and salinity are calculated using 
conservation equations, again utilizing a turbulent mixing hypothesis for closure. The 
equations are linked by a simplified equation of state. 

For the purpose of computational efficiency several techniques are used. High 
speed, external gravity waves are eliminated by the "rigid-lid" approximation, and a 
Laplacian equation is solved for the non-divergent, vertically averaged flow. The 
timestep limitation, i.e., the half pendulum day constraint associated with inertia-gravity 


waves, is overcome by a semi-implicit treatment of the Coriolis term. 


Considerable improvement was made to the structure of the FORTRAN code of 
this model by Semtner (1974) who, at the same time, added various features to the 
mathematical formulation, chief of which was the use of "hole relaxation" (Takano, 
1974) in the solution of the external mode for a model with islands. This version of the 
model has been adopted by many investigators and has seen considerable use for a 
number of years in the ocean modelling work at GFDL. During this time, as vector 
processing machines became more demanding of suitable FORTRAN structure, 
significant changes have been made to the code for efficiency purposes. It has also been 
generalized in several ways, among which is the incorporation of variable grid spacing 
in the horizontal, and an arbitrary number of tracer prognostic variables. The relaxation 
code for the solution of the external mode has been redesigned, and a better technique 
for establishing the initial guess has reduced the scans-to-convergence considerably. 

Two objectives have been sought in designing the code for use outside of GFDL. 
First, it has been made universal to some degree, by the use of optional code lines. A 


separate, FORTRAN coded updating utility is provided to carry out these operations. 


B. THE MODEL DESCRIPTION 

The GFDL Modular Ocean Model (acronym MOM) version 1.0...Dec 1, 1990 1s 
a three dimensional primitive equation general ocean circulation model intended for use 
as a flexible tool for exploring ocean and coupled air-sea applications over a wide range 


of space and time scales. 


MOM has been written as a collaborative effort by Ron Pacanowski, Keith Dixon, 
and Anthony Rosati at the National Oceanographic and Atmospheric Administration’ s 
Geophysical Fluid Dynamics Laboratory in Princeton, New Jersey. It is the successor 
to the code written by Michael Cox, documented in the GFDL Ocean Tech Report #1, 
(1984). As was the case for the Cox model and the Semtner model (UCLA Dept. of 
Meteorology Tech. Report No. 9, 1974) that preceded it, MOM is a Fortran 
implementation of equations described by Kirk Bryan (1969). 

For my investigation, I use the Pacific Ocean as the region of interest in the MOM 
model. The region is defined by 160° of longitude from 120°E to 280°E and extending 
from 60°S to 60°N. A constant depth of 5700 m is assumed. The horizontal resolution 
is 5° in longitude and 5° in latitude. There are 15 levels in the vertical, with resolution 
varying from 30 m near the surface to 836 m near the bottom. 

The time manager clock parameters for setting time of model initial conditions are 
year) = starting year (set to be 0), monthO = starting month (set to be 1), and dayO = 
starting day (set to be 1). Using logical Julian calendar in MOM model, the time period 
of integrations is 60 years. The time step for both barotropic and baroclinic velocity is 
l hour. For density and tracers A temperature and salinity), the time step is | 
day. 

The MOM model uses an efficient grid system. Horizontally, tracer quantities are 
defined at the centers of "t" grid boxes and velocities are defined at the centers of "u,v" 
grid boxes. The centers of "u,v" grid boxes are located at the northeast corners of "t" 


grid boxes. The first "t" grid box is located in the southwest corner of the "t" grid. 


This grid system is replicated and stacked vertically one on top of another from the 
surface of the ocean downward. Vertically, tracers and velocities are defined at the 
centers of their corresponding boxes and are set at the same depths. 

The values for the mixing coefficients of MOM model were chosen as 10? cm?/s 
and 1 cm”/s for the horizontal and vertical viscosity, respectively, and 2 X 10* cm?/s and 
20 cm’/s for the horizontal and vertical diffusivity, respectively. The bottom drag 


coefficient was set to Q. 


C. BASIC EQUATIONS OF MOM MODEL 


The basic equations of the model are written here in continuous form. 


Let 
m = sec 
n= sind (6) 
f = 2Qsinọ 


where @ is latitude. 


An advective operator, 


P (1) = m= [> (wH), (7) 





(UH), + [va 


is adopted in which p is any scalar quantity, À is longitude and a is the radius of the 
earth. 


The equations of motion are then 


Ve ae ie) Siku = E + FY 
b 


where pọ is taken to be unity. 


The local pressure, p, is given by the hydrostatic relation, 
0 
p(z) = p* + | gp dz 


where p° is the pressure at the surface of the ocean. 


The continuity equation is 


The conservation equation, 


(8) 


(9) 


(10) 


(ET) 


(12) 


applies to any "tracer" type of quantity carried in the model. These include the active 


tracers, potential temperature and salinity (active in the sense that they appear in the 


equation of state), and any passive tracers such as Carbon 14 or Tritium. 


The equation of state is 


p = pO 


(13) 


where @ is potential temperature, S is salinity and the depth dependence arises from 
compression effects. In the present model, (13) is represented by a polynomial fit to the 
Knudsen formula (Bryan and Cox, 1972). 


Let 


NA = má + 1 Le (14) 
> 


The effects of turbulent mixing are 


FY = Ayu, + Aga (Vou a i min eo (15) 
FY = AwV.,, + Awa “(Vv + (1 = m°n*) vga (16) 
A È 
zi- ae + Ana VER (17) 
zZ 





where A is the mixing coefficient corresponding to M (momentum), T (tracer), V 
(vertical), and H (horizontal). 

In nature, vertical mixing is known to be a rather complex function of vertical 
stability. Since this process is still not well understood, we have adopted a simple, 
uniform mixing under statically stable conditions, and an infinite mixing under statically 


unstable conditions. This is done by specifying 6 to be 








1 cp < 0 
e GZ (18) 
0 OP 5 9 
oz? 
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At lateral walls, the boundary conditions are 
Ua eee (19) 


where the nm subscript indicates a local derivative with respect to the coordinate normal 


to the wall. At the surface, 


Da eam) emi) 
w= 0 


The "“rigid-lid" assumption of zero vertical motion at the surface filters out high speed 
external gravity waves which would otherwise seriously limit the length of the time step 
used in the numerical integration. The quantities 7, 7? are the zonal and meridional 
components of the specified surface stress, and y is a flux through the surface, of the 
particular tracer involved. 


At the bottom, 


T_=0 (21) 


w= -mua"*H, - va H, 
Combining (8) and (9) with (10), we derive 


NI O WAY 
V = ue ~ a De 
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where 


= =u) A mga `f paz’ E (23) 
= = p / V 

ve = D o a ea + E aa 

Let us define 
= es Ue Vv = yee (25) 

where 
E = Hf paz (26) 
-H 


Then, since p' is not a function of depth, 


O = u - u: Vo ve (27) 


Since all terms on the right of (23) and (24) are known, (27) may be solved for the 
internal modes of momentum. Under the rigid-lid boundary condition, the external mode 


of momentum may be represented by a volume transport stream function, y, 
= = (ae War v = m(aH) *v, (28) 


This is shown by integrating (11) vertically, substituting (28) and noting that the 
boundary conditions (20) and (21) on w are satisfied. A prognostic equation for y may 
be obtained by averaging (22) vertically, and eliminating terms in p* by applying the curl 


operator, 


12 


17 77 = sI 
a a a 





re 29 
Vey (2 (29) 


Substituting (28), 




















Mal ,| Lol adv, (E) (30) 
Ha? |, mHa? ly Mm) 
The boundary condition on y at lateral walls, corresponding to (19) is 
Va = Y, =0 (31) 


This condition is satisfied by setting Y constant over each unconnected land mass 
comprising the ocean boundary. In the case of an enclosed basin with no islands, y may 
be set to zero over the boundary forming land mass. If, in addition, islands are present, 
the associated constant for each island reflects the net flow around the island and must 
therefore be predicted by the governing equations. The method used is "hole relaxation" 
in which the line integral of the quantity Vp*, taken around the island, is required to 
vanish. Averaging (22) vertically, integrating around the coast of the island and setting 


the contribution due to p° to zero, the predictive equation, 


ea À Hee lan - affa $ (Zal (32) 


is obtained. Applying the Stokes theorem yields a more useful form, 


eel alje e- 


Note that (33) is simply an area integral of (30), taken over the island. 


dA (33) 
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D. PROGRAM FLOW OF MOM MODEL 


The FORTRAN code of the model consists of the main program OCEAN, and 


seven subroutines. 


OCEAN: 


STEP: 


CLINIC: 


TRACER: 


RELAX: 


STATE: 


Their functions are described below. 
Performs all operations. This is done only once at the beginning of 
each run of the model. The program calls STEP once per timestep, 
and attends to operations which must be done at the end of each mun. 
Called once per timestep by OCEAN, it initializes various quantities, 
bootstraps the row-by-row computation of prognostic variables, 
manages the I/O for the latter, and performs various analysis 
procedures on the progressing solution. 
Called once per row by STEP, it computes the internal mode 
component of the u and vy velocities as well as the vorticity driving 
function for use by RELAX later in determining the external mode. 
Called once per row by STEP, it computes temperature, salinity, and 
any tracers which are carried in the model. 
Called once at the end of each timestep by STEP, it takes the vorticity 
driving function computed in CLINIC and, using sequential 
overralaxation, solves the Laplacian equation for the external mode of 
velocity in terms of a mass transport stream function. 
Called by CLINIC and TRACER once per row, and STEP in the 
bootstrap procedure, it computes normalized densities by using a third 


order polynomial fit to the Knudsen formula. 
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MATRIX: Called by STEP on specified timesteps, it is a general 2-dimensional 

array printing routine. 

ODAM: (Ocean Direct Access Manager) 

Used primarily by STEP, it is a set of routines which is responsible 
for handling the transfer of data between memory and disc (virtual 
disc residing in memory in the core contained mode). 

In summary, for one timestep, OCEAN calls STEP, which establishes the proper 
data in memory and calls CLINIC and TRACER row-by-row from south to north through 
the basin. Upon completing the final row, STEP calls RELAX and returns control to 
OCEAN, which may call STEP for another timestep. 

Two additional subroutines, FINDEX and FILTER are used when Fourier filtering 
is needed to overcome the timestep limitation arising from convergence of meridians at 


high latitudes. 


E. DISC VO SYSTEM 

Next to the equations themselves, the component of the code which accounts for 
the greatest complexity is the I/O system. The purpose of this system 1s twofold. First, 
having a complete record of all prognostic variables on permanent disc at the end of each 
timestep allows restarting an experiment from an earlier run, or from a machine 
malfunction. Secondly, during a run, it is generally impossible to fit the three 
dimensional arrays in memory entirely. The I/O system is designed to feed data to and 


from memory in a row-by-row manner. Data for one row, including all east-west and 
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vertical grid points is termed a "slab". At any one time, only the slabs necessary for the 
computation are present. While the computation for one row proceeds, the I/O system 
is feeding the slab just computed, back to disc and fetching the slab which will be needed 
to compute the next row. If the disc transfers are done fast enough, no wait for data will 
be needed at the beginning of each row and the system is said to be completely 
"buffered". 

Three disc units are needed for the process when using centered differencing in 
time; one for the N timestep data (read), one for the N-1 timestep data (read), and one 
for the newly computed N+1 timestep data (written). Since the N level on one timestep 
becomes the N-1 level for the next, it is convenient to permute the disc units to minimize 
data transfers. Thus, on timestep |, N-1 and N data are read from units 13 and 14, and 
the N+1 data is written to 15. On timestep 2, units 14 and 15 are read and the N+1 
data is written to 13. On timestep 3, 15 and 13 are read with the N+1 data going to 14, 
CLC. 

Possibly the most abstruse feature of the model is the manner in which the "slab 
incidental data” is handled. Just as the slab system reduces the row dimension of the 3- 
dimensional prognostic variables to 3, it can also be used to reduce the row dimension 
of 2-dimensional variables which would otherwise add considerably to the memory 
requirements. Also, if these variables are constant in time, there is no need to keep 
multiple time records of them. Consider two arrays, A and B, for which there is data, 
invariant in time, at each row and column horizontally across the grid. They may be 


carried in the model as "slab incidental data", thereby reducing their row dimension to 
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3. Furthermore, let A reside in the slab corresponding to the N-1 time level of the 
primary slab data, and B reside in the N level slab. The diagram below illustrates the 
6 timestep cycle which the permuting disc units execute, with N-1, N and N+1 denoting 


the primary slab data, and A, B denoting slab incidental data. 


TIMESTEP: l 2 3 > 5 6 7 

UNIT 

13 N-1,A N+1,B NB N-1,B N+LA NA NIA 
14 NB N1B N+1A NA N1A N+I1B NB 
15 N+1A NA N1A N+1B NB  N-1B N+lL,A 


Note that on even timesteps, the slab incidental data enters memory in the wrong 
slabs and must be switched between slabs N and N-1. Also, storing of A and B into the 
N+1, slab must be alternated on successive timesteps. In the base code, FKMU and 
WSY are "A" type arrays, and FKMT and WSX are "B" type arrays. 

Disc unit 12 is used primarily for the storage of 2-dimensional, horizontal fields. 
It is divided into 7+ blocks in the following manner: 

Blocks 1-3: — permuting blocks for the stream function at N-1, N, and N+1 time 

levels 

Block 4: reciprocal of total depth 


Blocks 5-6:  permuting blocks for former relaxation solutions 
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Block 7+: start and end indices (additional blocks are added as necessary 
when filtering indices fill block 7) 

Finally, unit 11 contains the timestep counter, total elapsed time, and the area and 
volume of the basin. 

Data is fed to and from these units by means of the entry points of subroutine 
ODAM. It, in turn, must use a direct access I/O package provided locally. If 
FORTRAN direct access I/O is available, with a facility to buffer the operations (such 
as a FIND statement), the QDAM calls in ODAM of the base code may be replaced by 
the appropriate FORTRAN statements. Otherwise, a specially written set of I/O utilities 
such as QDAM must be supplied. If, instead, the core-contained option is invoked, no 


such package in needed. 
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III. COMPATIBILITY BETWEEN SURFACE WIND STRESS AND 
BUOYANCY FORCING 


Wind stress and buoyancy forcing are independently designated in most ocean 
models, e.g., no constraint has been considered between the surface wind stress and 
surface air temperature (or net surface heat flux). The atmosphere actually is a unified 
physical system. For atmospheric motion on the scale of the ocean basin considered 
here, the wind and temperature fields are related to each other through the thermal wind 
relation. Only one of the two (either wind or air temperature) can be prescribed in the 
atmosphere-driven ocean models. Thus, treating the 2 fields as totally independent will 
produce physically unrealistic atmospheric forcing in the ocean models. This gives rise 
to a distortion of model results since the solutions largely depend upon the atmospheric 
forcing. Therefore, it is important to study the compatibility between the wind and 
buoyancy forcing and the effect of incompatibility on the model results before running 
any ocean numerical model. 

| In many ocean models (e.g., Davis, 1983, Robinson et al. 1977, etc.) the SSAT 
is assumed zonal symmetric. Therefore, if the surface wind (u,”, v,) and SSAT fields 
used in ocean models are compatible then they can be denoted by combinations of 


sinusoidal functions (Chu, 1992), 1.e., 


AO AS Kray (34) 


n 
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Vee a Vi, Sill ery (35) 
n 


O.. = (AT) Soe cos k,ny EVR) 


i 


where the subscript ’*’ indicates the dimensional variables, k, (n = 1, 2,...) are not 


necessary integers. It should be noted that (u,", v,) is the total surface wind, i.e., 


Vi ae ano (37) 


x7 S* x7 


Since the flow variables are independent of x (zonally symmetric), the geostrophic 


flow should be purely zonal, therefore, at the APBL top the geostrophic wind would be 
Unea O ae, = 20 (38) 
The scale of thermally forced surface winds (Chu, 1989) is 


-A JONT 


CON (39) 
AS 


where 6 = (v/)”, is the Ekman depth, v is the atmospheric vertical eddy viscosity, Q 
is the angular velocity of the earth’s rotation, and AT/L, is the characteristic SSAT 
gradient. If 6 = 1 km (Holton, 1972), AT = 14°k, L, = 2000 km (Robinson et al, 
1977), and 0, = 288°k, the scale for the thermally driven surface wind (U;) approaches 
3 m/s, and the scale of the geostrophic wind at the APBL top (U) can reasonably be 
assumed to be 10 m/s. Therefore, if we completely neglect this thermally driven surface 


wind in ocean models, we will make almost 30% error in the surface wind forcing. 


20 


Based on the K-closure APBL model proposed by Kuo (1973), Chu (1992) found 


these wind and temperature fields are: 


uly,z) = Y [-b, exp(-y2Y 2) 


3 


a À 
+ 2f Ww o O 
m E oa A g (40) 


E y2Y fyb, 2 
O(y,z) = )) [~~ “exp (-v2y 2z) 


+ 


a, 
2 ma y se ay | 4,32) ] cos k,ny 


LAR, -zy (41) 


v(y, Z) ee E AnjAnjexp ( (A, 32) sin K,Ty (42) 


evaluated at z=Q, 1.e., 


2 a À 
O DE Cer PCa AA isn (43) 
J = 2 An = 
3 
on. + y a (44) 
al 
YA ano 
0. = MI 2k,nRi y _=— (45) 
KT J= ae 2Y 


where y is the nondimensional decay rate (y = 0.01, denoting 100 day decay time scale) 
for Rayleigh friction and Newtonian cooling, fọ = sinó, is the nondimensional Coriolis 


parameter, taken as locally constant except very close to the equator, @ is latitude, and 


2A 


Ri = Ri(8/H,)? (46) 


is the generalized Richardson number. H, = 10 km, is the thickness, and 


Ri =a (47) 


is the Richardson number. Here, N is the Brunt-Vaisala frequency. A,{j = 1, 2, 3) are 
the eigenvalues. a,,, a,,, and a,; are integral constants. y is the streamfunction. 
The boundary condition in the vertical are listed as follow. The dependent 


variables should remain finite as z > œ, i.e., 


y lul. J8)) < (48) 





du Ov 08 
Y|,.,=0, HU, + U = a VEZ TA 8 - eae = 0: (49) 
where M, M, are given constants, and 
u = U/Up = 3.33 (50) 


the parameter M = 0.5, which is a measure of the effective depth of the constant stress 
sub-layer (Kuo, 1973). 


The surface boundary condition (48) and (49) leads to (assuming M, = 0) 


$ ay =O (51) 
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Solving the linear algebraic equation (44), (51), and (53) with respect to a,,, a,,, and a,;, 


we obtain 
l l J 
a m a E (54) 
ni 1, sn n2 m sn na JA sn 
where 
0 1 JE 
7 = i NES À n3 (55) 


ni 


1/M Año Ms 


T 0 J 
L = À n1 1 À n3 (56) 
Aar IIM Na 
1 T 0 
ae = À n1 À n2 1 (57) 
Anz Ang 1/M 
ile 1 Jl 
Ll = Ani Anz À n3 (58) 


Ana An2 Ans 
Elimination of 6, from (43), (45), and (52) leads to 
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Un Sue lee eee (595 


sn fh sr 








Ven ~ -Rna E (60) 
where 

kT 
KE AI = (61) 

CE 

K TM 
Rin = ees F (62) 
0 

L 

Ro: HF Ha, = = (63) 
Nn n 
where 
2Mk£1?Ri z ee 
I, = (——~——- - 4Mf,,7) ` rr (64) 
fo JUL 1,(44, a 
2k*n?Ri 3 za 
Si A + /2y M) - 4f,My y —— 
£542Y j=1 lAn 7 2V) 

+ 2f — E 

02, ian - 2y) (65) 


The solutions (59), and (60) show that the surface winds (u,,, v,,) are driven 


sn? 


mechanically by the geostrophic winds, u,,, and thermally by the surface thermal 


conditions, -R,,0.,, and -R,,6,,. 


The coefficient in (59) and (60) R,,, R»,, H,,, and H,, for different k, can be 


computed by using the dependence of eigenvalues A,,, A,2, and A,; On k, for 
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eae Ot (66) 


Non-zero values for these coefficients indicate that only two can be prescribed 
among the four variables u,,, V,,, 9,,, and uc, , the other two should be derived from (59) 
and (60). Thus, (59) and (60) can be treated as a compatibility condition among these 
four variables. Obviously, it is over-specified if we prescribe both surface wind stress 
and buoyancy forcing in ocean numerical models. 

We can use (59) and (60) to verify the consistency between buoyancy forcing (3) 
and surface wind (4) which is used in the Davey (1983) model. Here, 7,(y) = O, the 


Fourier cosine series of 7, is 


Bae = (il, = POD o7 COS NIY (67) 


The surface winds (u,, v,) are set to zero, therefore, the components of their Fourier sine 


series should also be zero, 1.e., 


eel, Via el (68) 
Non-zero values of 
p, = ee E UPR) (69) 
A 
2n 


shows that the linear algebraic equations (59) and (60) about @,, and u,, only have zero 
solutions, i.e., 6, = 0, ue, = 0, when u,, = 0, v,, = 0. Therefore, the thermal forcing 


’ sn 


(3) is not compatible with the surface wind forcing (4). 
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Purely zonal surface winds are commonly used in ocean models, such as the 


$2 282 


surface boundary conditions (5). Here, three variables u,, v,, and J, are prescribed. We 


expand these surface variables into Fourier sine and cosine series 





T 
(Ugr ve = : e MD View sin lea, Vian = O (70) 
P2Cp n 
ie = (AT eon COS nny (719 
n 
UH ASIA (72) 


n 


The Fourier components tsn, Vens Osn, and lg, should satisfy the two linear algebraic 


equations (59) and (60): 
Usn = SR O ja H nUgn (73) 


OS Sh eee ja Hz nUgn (74) 


which shows that purely zonal surface winds appear only when the surface thermal 


forcing is balanced by the geostrophic forcing 





0. = R Ugn (75) 


Elimination of u¿, from these two equations leads to a relationship between u,, and 


9... (Chu, 1992): 


Usn = BS (76) 
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which is the compatibility condition for zonally symmetric zonal wind and SSAT fields. 
Independently prescribing surface wind and SSAT fields, such as in (5), violates this 


condition and therefore, is inconsistent. 
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IV. ERRORS IN NUMERICAL INTEGRATION OF THE MOM MODEL 
CAUSED BY INCOMPATIBLE FORCING 
A. SURFACE BOUNDARY CONDITIONS 

The MOM model and its parameters have been discussed in Chapter II. Two kinds 
(compatible and incompatible) of surface boundary conditions (SST, salinity, wind stress 
in x direction, and wind stress in y direction) were used in MOM model to find the 
differences between the results of the two conditions. 

The incompatible surface boundary condition (SST, salinity, wind stress in x 
direction, and wind stress in y direction) which is used in the general circulation model 
by Marotzke and Willebrand (1991) is shown in figure la. The zonal (eastward) wind 
stress (solid curve) and the apparent atmospheric temperature in degrees celsius (dashed 
curve) are functions of latitude and symmetric about the equator, so they are only shown 
for one hemisphere. 

The sea surface temperature (SST) in °C is symmetric about the equator, and 
follows a cosine law in latitude with a difference of 27°C between high and low latitudes 
(Fig. 1b). The salinity is analogous to the sea surface temperature and symmetric about 
the equator. The salinity also follows a cosine law in latitude with amplitude 2.5 psu 


(Fig. Lc). 
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Many independently and simple prescribed wind and buoyancy forcing in ocean 
models have a similar form as used by Marotzke and Willebrand (1991) in a steady-state 
wind and thermally forced ocean circulation model (Fig. 1, 2 and 3). 

The latitudinal surface atmospheric temperature gradient generates extra surface 
winds (Up Vr). The scale of this thermally driven surface winds is given by (39). 
Corresponding to SSAT field described in figure 1b, the thermally driven surface winds 
can be predicted by a K-closure Atmospheric Planetary Boundary Layer (APBL) model 


(Chu, 1989). From wind and temperature fields (40), (41) and (42), we obtain 


T e 
7 Ur COSA -Ssina | ji (77) 
G,|sine cose sere" 





Vr 








where a is the deflection angle (angle between V; and V@, | ,-.). The wind field with 
the thermally driven component (i.e., u+u;, v7) is illustrated by solid curves and without 
the thermally driven component is illustrated by dashed curves in figure 2 and 3. 

From figure 2 and 3, we compute the Ekman current change due to compatible 
forcing (Fig. 4). The result in a stronger westward and poleward current between about 
5° and 30° latitude in both hemisphere. We also compute the wind stress curl for the 
wind field with and without thermally driven component. 


Wind stress curl for forcing without thermally driven component: 


R-curi¢ = Sx - Ss (78) 
xX Ja 
7, 18 independent of x, so 
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E a 
K- Cur! c Aa (79) 
dy 


Wind stress curl for forcing with thermally driven component: 


— Ot za 
K Cor A e (80) 
Ox oy 
z, is independent with x, so 
A OT 
E Curl AA (81) 
dy 


The difference between the wind stress curl for wind field with and without thermally 


driven component is not negligible. 


B. ERRORS CAUSED BY INCOMPATIBLE WIND AND BUOYANCY 
FORCING 
We consider a rectangular ocean of uniform depth in the MOM model. The 
thickness (m) and depth (m) in the vertical of the 15 level MOM model is specified as 


below: 


76.15 
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1. The total northward transport of heat 











The total northward transport of heat, which was calculated by compatible 
surface wind stress and buoyancy forcing, after 3, 6, 9, and 12 years is shown in figure 
5. The difference in the total northward transport of heat between compatible and 
incompatible surface wind stress and buoyancy forcing after 3, 6, 9, and 12 years is 
provided in figure 6. In all plots the difference is computed from the compatible forcing 
run minus the incompatible forcing run. 

The total northward transport of heat, which was calculated by compatible 
surface wind stress and buoyancy forcing, after 51, 54, 57, and 60 years is shown in 
figure 7. The difference in the total northward transport of heat between compatible and 
incompatible surface wind stress and buoyancy forcing after 51, 54, 57, and 60 years is 


provided in figure 8. 
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When the northward transport of heat by compatibility consideration is large 
in figure 5, the difference in the northward transport of heat, shown in figure 6, due to 
the incompatibility forcing, is also large. The error reaches to 21% after 51 year’s 
integration. The difference in the northward transport of heat between the compatibility 
and incompatibility consideration increases with the integration time. 

In northern (southern) hemisphere, an additional northeasterly (southeastly) 
wind stress exists in the compatible forcing case (Fig. 2, 3). Therefore, an extra 
northwest (southwest) Ekman transport should appear when we use the compatible 
forcing (Fig. 4). 

2. The total northward transport of salt 

The total northward transport of salt which was calculated by compatible 
surface wind stress and buoyancy forcing, after 3, 6, 9, and 12 years is shown in figure 
9. The difference in the total northward transport of salt between compatible and 
incompatible surface wind stress and buoyancy forcing after 3, 6, 9, and 12 years is 
provided in figure 10. 

The total northward transport of salt which was calculated by compatible 
surface wind stress and buoyancy forcing, after 51, 54, 57, and 60 years is shown in 
figure 11. The difference in the total northward transport of salt between compatible and 
incompatible surface wind stress and buoyancy forcing after 51, 54, 57, and 60 years is 
provided in figure 12. 

When northward transport of salt by compatibility consideration is large in 


figure 9, the difference in the northward transport of salt in figure 10, due to the 
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incompatibility forcing, is also large. The error reaches to 16% after 51 year’s 
integration. The difference in the northward transport of salt between the compatibility 


and incompatibility consideration increases with the integration time. 


3. The meridional mass transport 

The meridional mass transport, which was calculated by compatible surface 
wind stress and buoyancy forcing, after 30 years is shown in figure 13. The difference 
in the meridional mass transport between compatible and incompatible surface wind stress 
and buoyancy forcing after 30 years is provided in figure 14. 

The meridional mass transport, which was calculated by compatible surface 
wind stress and buoyancy forcing, after 60 years is shown in figure 15. The difference 
in the meridional mass transport between compatible and incompatible surface wind stress 
and buoyancy forcing after 60 years is provided in figure 16. 

From above, the latitudinal distribution of meridional mass transport after 60 
years integration indicates that the compatibility consideration meridional mass transport 
is larger than that of incompatibility consideration in all latitudes of the northern 
hemisphere and smaller in all latitudes of southern hemisphere. The difference between 
the two consideration (compatibility and incompatibility) increases with the latitude. 

In figure 16, the difference of mass transport is strong between 25°S and 
25°N due to the extra Ekman transport when the compatible forcing is used (Fig. 4). 
The error is more then 5 Sverdrup in northward mass transport. In deeper sea water it 


is about 5 Sverdrup less. 
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4. Meridional section profile of 0 

The meridional section profile of O at 160%E after 60 years is shown in figure 
17 and the difference is provided in figure 18. The meridional section profile of 0 at 
160% W after 60 years is shown in figure 19 and the difference is provided in figure 20. 
The meridional section profile of O at 120%W after 60 years is shown in figure 21 and 
the difference is provided in figure 22. 

From the difference of meridional section profile of @ in figure 18, 20, and 
22, we can see that the potential temperature @ 1s always lower in the compatible forcing 
than in incompatible forcing near the equator. The errors caused by incompatible forcing 
in upper is larger than that in the deeper sea water. 

The mass transport is northward in north hemisphere and southward in south 
hemisphere. Thus divergence occurs near the equator. So the cold deeper sea water 


moves upward toward the surface to make the temperature lower near the equator. 


5. Meridional section profile of v velocity 
The meridional section profile of v velocity at 160°E after 60 years is shown 
in figure 25, and the difference is provided in figure 26. The meridional section profile 
of v velocity at 160°W after 60 years is shown in figure 27, and the difference is 
provided in figure 28. The meridional section profile of v velocity at 120°W after 60 
years is shown in figure 29, and the difference is provided in figure 30. 
From the meridional section profiles of v velocity in figure 25, 27, and 29, 


the v velocity is positive from the equator to about 25°N where the velocity changes 
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direction (opposite sign) toward high latitude. The v velocity in the southern-hemisphere 
is compared with the relation in the northern-hemisphere. 

From figure 14 and 16, the latitudinal distribution of meridional mass 
transport after 60 years integration indicates that the compatibility consideration 
meridional mass transport is larger than that of incompatibility consideration in all 
latitudes of the northern hemisphere and smaller in all latitudes of southern hemisphere. 
This is consistent with figure 4, which shows an extra northwest (southwest) Ekman 
transport in northern (southern) hemisphere when the compatible forcing is used. 

Therefore, we see the difference of meridional section profile of v velocity in 
figure 24, 26, and 28, the v velocity, which is due to compatible forcing correction, is 
always stronger in the northern-hemisphere than that due to incompatible forcing 
influence. In the southern-hemisphere the resulting v velocity which is integrated by 
compatible forcing is found to be weak. The errors by incompatible forcing is most 


obvious in the upper sea water. The errors is nearly 25%. 


6. Meridional section profile of w velocity 
The meridional section profile of w velocity at 160°E after 60 years is shown 
in figure 29, and the difference is provided in figure 30. The meridional section profile 
of w velocity at 160°W after 60 years is shown in figure 31, and the difference is 
provided in figure 32. The meridional section profile of w velocity at 120°W after 60 
years is shown in figure 33, and the difference is provided in figure 34. 
From the meridional section profiles of w velocity in figure 29, 31, and 33, 


the w velocity has large variation in longitude (i.e., section to section) indicating the 
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effect of transients (eddies). The w velocity profile of deep water is more complex than 
that of upper sea water. From the difference of meridional section profile of w velocity 
in figure 30, 32, and 34, the w velocity, which is due to incompatible forcing influences, 


has more errors near the equator than that of high latitude. The errors is about 16%. 
7. Difference at some random points 


a. 0 at 160°E, O°N on level 2 

The difference in @ for depth of 76.15 meters (level 2) at 160°E, O°N 
is shown in figure 35. The error for @ in the first three years is about 6% due to 
incompatible forcing influences. This error increases to 7% after about 10 years 
integration. 

An extra northwest (southwest) Ekman transport appears in the northern 
(southern) hemisphere when the compatible forcing is used. Thus an additional 
horizontal divergence occurs near the equator. Causing an extra amount of cold deeper 
sea water moves upward toward the surface making the temperature lower near the 
equator in the compatible case. Therefore, at depth of 76.15 meters (level 2) at 160°E, 
O°N the @ by the compatible forcing in the whole integration time of 60 years is always 


lower than that due to incompatible forcing influences. 


b. V velocity at 160°E, 5°S on level 2 
The difference in v velocity for depth of 76.15 meters (level 2) at 160°E, 


5°S is provided in figure 36. The error in v velocity in the first three years is about 
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26% due to incompatible forcing influences. This error increases to about 28% after 
about 8 years integration. 

From the meridional section profiles of v velocity in figure 23, 25, and 
27, the v velocity is positive from the equator to about 25°N where the velocity changes 
direction (opposite sign) toward high latitude. The v velocity in the southern-hemisphere 
is compared with the relation in the northern-hemisphere. 

From figure 14 and 16, the latitudinal distribution of meridional mass 
transport after 60 years integration indicates that the compatibility consideration 
meridional mass transport is larger than that of incompatibility consideration in all 
latitudes of the northern hemisphere and smaller in all latitudes of southern hemisphere. 
This is consistent with figure 4, which shows an extra northwest (southwest) Ekman 
transport in northern (southern) hemisphere when the compatible forcing is used. 

Therefore, at depth of 76.15 meters (level 2) at 160°E, 5°S the v 
velocity by the compatible forcing in the whole integration time of 60 years is always 


weaker than that due to incompatible forcing influences. 


c. V velocity at 160°E, 10°S on level 2 
The difference in v velocity for depth of 76.15 meters (level 2 at 160°E, 
10°S can be seen in figure 37. The error in v velocity in the first three years is about 
23% which is due to incompatible forcing influences. This increases the error to 28% 
after about 8 years integration. At depth of 76.15 meters (level 2) at 160°E, 10°S the 
v velocity by the compatible forcing in the whole integration time of 60 years is always 


weaker than that due to incompatible forcing influences. 
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d. W velocity at 160°E, 10°N on level 5 
The difference in w velocity for depth of 385.64 meters (level 5) at 
160°E, 10°N is illustrated in figure 38. The error of w velocity in the first three years 
is about 16% due to incompatible forcing influences. This error decreases to 10% after 
about 8 years integration. At depth of 385.64 meters (level 5) at 160°E, 10°N the w 
velocity by the compatible forcing is always stronger than that due to incompatible 


forcing influences after 6 years integration. 
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Figure la. The surface forcing fields which were used by Marotzke and 
Willebrand (1991). 


-60 -50 40 -30 -20 -10 0 10 20 30 40 50 60 
LATITUDE 





Figure 1b. Zonal mean values of observed sea surface temperature (*C), used by 
- Marotzke and Willebrand (1991). 
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Figure lc. Zonal mean values of observed salinity (psu), used by Marotzke and 
Willebrand (1991). 
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Figure 2. Surface zonal wind stresses (dynes/cm’) as functions of latitude. 
(a) without thermally driven surface winds (dash curve). (b) with 
thermally driven surface winds (solid curve). 
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Figure 3. Surface meridional wind stresses (dynes/cm’) as functions of 
latitude. (a) without thermally driven surface winds (dashed 
curve). (b) with thermally driven surface winds (dashed curve). 
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Figure 4. The Ekman current change near 25°N and 25°S due to 
compatible forcing. The solid arrows show the change in the 
wind stress (Figs. 2 and 3) and the double arrows show the 
resulting change in the surface Ekman current. 
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DASH: After 6 years 
DOT: After 9 years 
DASH-DOT: After 12 years 





Figure 5. Total northward transport of heat after 3, 6, 9, 12 years 
integration. (x 105 Watts). 
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Figure 6. Difference of total northward transport of heat after 3, 6, 9, 12 
years integration. (x 10° Watts). 
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SOLID: After 51 years 
DASH: After 54 years 
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Figure 7. Total northward transport of heat after 51, 54, 57, 60 years 
integration. (X 10% Watts). 
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Figure 8. Difference of total northward transport of heat after 51, 54, 57, 
60 years integration (X 10% Watts). 
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SOLID: After 3 years 
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Figure 9. Total northward transport of salt after 3, 6, 9, 12 years 
integration. (x 10'° cm?/sec). 
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Figure 10. Difference of total northward transport of salt after 3, 6, 9, 12 years 
integration. (X 10!° cm?/sec). 
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Figure 11. Total northward transport of salt after 51, 54, 57, 60 
integration. (X 10' cm?/sec). 
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Figure 12. Difference of total northward transport of salt after 51, 54, 57, 60 


years integration. (x 10'° cm?/sec). 
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Meridional mass transport after 30 years integration. (in Sverdrup). 


Figure 13. 


Difference of meridional mass transport after 30 years integration. (in 


Sverdrup). 





Figure 14. 
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Meridional mass transport after 60 years integration. (in Sverdrup). 


Figure 15. 





Difference of meridional mass transport after 60 years integration. (in 


Sverdrup). 


Figure 16. 
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Figure 17. Distribution of temperature (°C) at 160°E meridional section after 60 
years integration. 





Figure 18. Distribution of temperature difference (°C) at 160°E meridional 
section after 60 years integration between two different kinds of 
surface forcing. 
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Figure 19. Distribution of temperature (°C) at 160°W meridional section after 60 
years integration. 





Figure 20. Distribution of temperature difference (°C) at 160°W meridional 
section after 60 years integration between two different kinds of 
surface forcing. 
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Figure 21. Distribution of temperature (°C) at 120°W meridional section after 60 
years integration. 





Figure 22. Distribution: of temperature difference (°C) at 120°W meridional 
section after 60 years integration between two different kinds of 
surface forcing. 
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E meridional section after 60 years 


Distribution of v (cm/sec) at 160° 


integration. 


Figure 23. 





Distribution of v difference (cm/sec) at 160°E meridional section after 
60 years integration between two different kinds of surface forcing. 


Figure 24. 
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Figure 25. Distribution of v (cm/sec) at 160°W meridional section after 60 years 
integration. 





Figure 26. Distribution of v difference (cm/sec) at 160°W meridional section 
after 60 years integration between two different kinds of surface 
forcing. 
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Figure 27. Distribution of v (cm/sec) at 120°W meridional section after 60 years 
integration. 





Figure 28. Distribution of v difference (cm/sec) at 120°W meridional section 
after 60 years integration between two different kinds of surface 
forcing. 
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Figure 29. Distribution of w (X 10° cm/sec) at 160°E meridional section after 60 
years integration. 





Figure 30. Distribution of w difference (x 10° cm/sec) at 160°E meridional 
section after 60 years integration between two different kinds of 
surface forcing. 
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Figure 31. Distribution of w (X 10° cm/sec) at 160°W meridional section after 
60 years integration. 
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Figure 32. Distribution of w difference (Xx 10? cm/sec) at 160W meridional 
section after 60 years integration between two different kinds of 
surface forcing. 
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Figure 33. Distribution of w (X 10° cm/sec) at 120°W meridional section after 
60 years integration. 
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Figure 34. Distribution of w difference (X 10° cm/sec) at 120°W meridional 
section after 60 years integration between two different kinds of 
surface forcing. 
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Figure 35. 
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Difference in 6 (°C) for depth of 76.15 meters (level 2) at 160%E, 
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Figure 37. Difference in v (cm/sec) for depth of 76.15 meters (level 2) at 160%E, 
1095. 





Figure 38. Difference in w (X 10? m/day) for depth of 385.64 meters (level 5) 
at 160°E, 10°N. 
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V. CONCLUSIONS 

The surface wind and surface thermal conditions (e.g., SSAT, net surface heat flux) 
are important components of the atmospheric system and are related each other. 
However, the atmospheric parameters representing mechanical and buoyancy forcing, 
such as surface wind and SST fields (or the net surface heat flux) are often prescribed 
independently in current ocean models, especially idealized or "process" models. The 
practice of prescribing wind stress and buoyancy forcing independently in ocean 
numerical models is an unrealistic assumption. 

A large difference between the ocean circulation forced by compatible and 
incompatible surface forcing (1.e., surface wind stress and buoyancy forcing) has been 
found by comparing the results from the GFDL MOM model. In particular, the 
compatible forcing derived from the (incompatible) forcings field of Marotzke and 
Willebrand (1991), produces 15-20% stronger meridional circulation cells in the low 
latitudes of both hemispheres of the MOM model. This difference can be attributed to 
the incompatible SSAT (or surface heat flux) and surface wind fields. Since most ocean 
numerical models are forced models, a compatible forcing boundary condition is an 


importance consideration before running any numerical ocean model. 
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